TP 2 Clasificacion. Tomas Palazzo y Axel Fridman
Cargamos las librerias
library("ggplot2")
library("GGally")
df = read.csv("lluviaAus.csv")
Seteo de semilla random
set.seed(329)
Limpieza y pre procesamiento de datos:
df$RainToday = as.factor(df$RainToday) # paso variables categoricas como factor
df$RainTomorrow = as.factor(df$RainTomorrow)
df$X <- NULL # borro columna X ya que sospechamos que no representa nada sino que es algun tipo de "id" que quedo grabado en el dataframe y no tiene influencia en la observacion.
Chequeamos que cada columna sea del tipo correcto
str(df)
'data.frame': 1000 obs. of 12 variables:
$ MinTemp : num 17.8 9 7.8 6.5 9 17.4 18.6 18.9 16.4 8.7 ...
$ MaxTemp : num 25.2 16 12.2 17.5 22.6 33.4 32.6 35.5 22.9 22.1 ...
$ Rainfall : num 0 0.8 1.8 0 0 0 0 0 0 0 ...
$ Evaporation : num 4 1.6 0.6 2 2.8 6.8 9 15.2 3.6 3 ...
$ Sunshine : num 6.4 7.4 5.4 9.7 9.5 10.5 11.8 12.2 10.4 8 ...
$ WindSpeed3pm: int 13 26 22 7 11 20 15 31 13 15 ...
$ Humidity3pm : int 66 53 84 47 42 18 38 13 57 39 ...
$ Pressure3pm : num 1013 1013 1026 1025 1016 ...
$ Cloud3pm : int 7 7 6 2 5 1 1 1 2 5 ...
$ Temp3pm : num 24.4 14.8 10.5 17.2 21.6 32.1 29.8 34.9 21.8 21.5 ...
$ RainToday : Factor w/ 2 levels "No","Yes": 1 1 2 1 1 1 1 1 1 1 ...
$ RainTomorrow: Factor w/ 2 levels "No","Yes": 1 1 2 1 2 1 1 1 1 1 ...
Divido dataset en train y test
sample <- sample(c(TRUE, FALSE), nrow(df), replace=TRUE, prob=c(0.8,0.2))
dftrain <- data.frame(df[sample, ])
dftest <- data.frame(df[!sample, ])
Analisis exploratorio de datos y visualizaciones:
g = ggpairs(df, progress = FALSE, bins=10)+theme_bw()
g

Lo que observamos es que hay variables que estan muy correlacionadas
linealmente, como la de temperatura maxima y temperatura a las 3pm. Y
otras que parecen ser muy poco relacionadas como la humedad y la
presion.
print(paste("Correlacion humedad y presion (baja) " , cor(df$Humidity3pm, df$Pressure3pm) ))
[1] "Correlacion humedad y presion (baja) 0.0456887832592021"
print(paste("Correlacion maxTemp y temp a 3pm (muy alta)" , cor(df$Temp3pm, df$MaxTemp) ))
[1] "Correlacion maxTemp y temp a 3pm (muy alta) 0.980580835987693"
Tambien notamos que aproximadamente 4/5 de las observaciones no
llueve, tanto ese mismo dia como el siguiente.
table(df$RainToday)
No Yes
802 198
table(df$RainTomorrow)
No Yes
794 206
Pero no son para nada independientes si llueve hoy y si llueve
mañana. Ya que, si solo tuviera de informacion estas 2 columnas, dado
que llovio hoy la nueva probabilidad de que llueva mañana es
aproximadamente 65/(97+65) =aprox 40% mucho mas que el 20% naive.
dfLluvias = df[c("RainToday", "RainTomorrow")]
table(dfLluvias)
RainTomorrow
RainToday No Yes
No 674 128
Yes 120 78
Ejercicio 2
ggplot(df, aes(x=Sunshine, y=Humidity3pm, color=RainTomorrow)) +
geom_point()

Observamos cosas muy relevantes, los dias que va a llover mañana,
tienen mayor humedad y menos sol, y los dias que no llovera mañana
tienen todos mucho mas sol y la humedad en promedio es mas baja. A su
vez hay varios dias, en su mayoria dias que llovera mañana, cuyo nivel
de sol es 0, lo cual genera esa columna en el lado izquierdo. Tambien
vemos que si bien hay muchos dias que tienen mucho sol y poca humedad
(los dias que no llovera mañana), no vemos casi ninguna observacion con
poco sol y poca humedad, lo cual nos podria hablar de cierta relacion
humedad - sol. Nosotros pensamos que la humedad tiene mayor capacidad
predictiva (si se tomara una sola y no en conjunto ambas), ya que la
variable sunshine esta mucho mas dispersa para los dias que No llovio al
dia siguiente, lo analizaremos en dos graficos de densidad.
ps<-ggplot(df, aes(x=Sunshine, fill=RainTomorrow)) +
geom_density(alpha=0.4) + labs(x= "Nivel de radiacion solar (Sunshine)",
subtitle="Grafico densidad radiacion solar") + geom_vline(xintercept=7.8, size=0.5, color="red")
ps

ph<-ggplot(df, aes(x=Humidity3pm, fill=RainTomorrow)) +
geom_density(alpha=0.4) + labs(x= "Nivel de humedad (Humidity3pm)",
subtitle="Grafico densidad humedad") + geom_vline(xintercept=62, size=0.5, color="red")
ph

Despues de ver estos 2 graficos notamos que no es facil dar un punto
de corte para diferenciar las 2 clases solamente tomando una variable.
Ya que si tomasemos aproximadamente 7.8 de punto de corte para la
radiacion solar o 62 para la humedad como punto de corte, de todas
formas tendrias bastante error ya que las clases se solapan mucho
mirandolas unidimensionalmente. Tomamos estos valores de referencia como
para decir que ningun corte es bueno, estos ni siquiera tienen en cuenta
la diferencia de proporcion de clases. En definitiva no nos casamos con
ninguna variable.
Ejercicio 3. Los boxplot no son buenos graficos. Pueden ocultar
demasiada informacion cuando hoy en dia tenemos la capacidad de
procesarla.
p2<-ggplot(df, aes(y=Humidity3pm, x=RainTomorrow, fill=RainTomorrow)) +
geom_boxplot() + labs(x= "Nivel de humedad (Humidity3pm)",
subtitle="Boxplots humedad segun lluvia mañana")
p2

p3<-ggplot(df, aes(y=Sunshine, x=RainTomorrow, fill=RainTomorrow)) +
geom_boxplot() + labs(x= "Nivel de radiacion solar (Sunshine)",
subtitle="Boxplots radiacion solar segun lluvia mañana")
p3

Como ven las medias difieren para ambos casos, pero esa informacion
ya la habiamos visto (ademas de muchas otras cosas que aca no) en los
density plots. No hay outliers / valores atipicos.
Ejercicio 4 Para hacer las ventanas moviles voy a primero transformar
el dataset en 1 y 0 a las categoricas, para poder luego tomar
promedios.
dftrain$RainToday = ifelse(dftrain$RainToday=="Yes",1,0)
dftrain$RainTomorrow = ifelse(dftrain$RainTomorrow=="Yes",1,0)
dftest$RainToday = ifelse(dftest$RainToday=="Yes",1,0)
dftest$RainTomorrow = ifelse(dftest$RainTomorrow=="Yes",1,0)
promediosMoviles = function(datosX, datosY, valor, h){
dfe = data.frame(datosX,datosY)
df2 = dfe[dfe$datosX >= valor-h & dfe$datosX <= valor+h ,]
while(nrow(df2)==0){
h = h*2 # Si no me agarra a nadie para promediar, aumentame la ventanita. No suele pasar, solo outliers o extremos
df2 = dfe[dfe$datosX >= valor-h & dfe$datosX <= valor+h ,]
}
if(mean(df2$datosY)>1/2 ){
return(1)
}
return(0)
}
promediosMoviles(dftrain$Sunshine, dftrain$RainTomorrow, 8, 0.1)
[1] 0
promediosMoviles(dftrain$Humidity3pm, dftrain$RainTomorrow, 85, 2)
[1] 1
Ejercicio 5 Nos creamos una funcion que nos genere todo el vector de
predicciones Yhat.
Vamos a evaluarlo con el metodo de validacion cruzada de LOO (dejar
uno afuera para entrenar y evaluarlo con ese).
leaveOneOut = function(datosX, datosY, h){
error = 0
for (i in (1: length(datosX))){
predichoI = promediosMoviles(datosX[-i], datosY[-i], datosX[i], h)
error = error + abs(predichoI - datosY[i])
}
return(error)
}
hPosibleHumedad = seq(1, 30, 0.5 )
hPosibleSunshine = seq(0.5, 10, 0.1 )
erroreshHumedad = c()
hum = dftrain$Humidity3pm
lluv = dftrain$RainTomorrow
for (i in (1: length(hPosibleHumedad))){
erroreshHumedad[i] = leaveOneOut(hum, lluv, hPosibleHumedad[i])
}
erroreshSunshine = c()
hum = dftrain$Humidity3pm
lluv = dftrain$RainTomorrow
sun = dftrain$Sunshine
for (i in (1: length(hPosibleSunshine))){
erroreshSunshine[i] = leaveOneOut(sun, lluv, hPosibleSunshine[i])
}
plot(hPosibleHumedad , erroreshHumedad, type = "l")

De aca vemos que la ventana optima para humedad es 1.
plot(hPosibleSunshine , erroreshSunshine, type = "l")

Mientras que la ventana optima para la radiacion solar es de 0.8
Ejercicio 6 ahora nos fijamos el error con el dataset de test que
separamos al principio para esta ventana
generarColumnaPrediccionesPromediosMoviles = function(datosX, datosY , h){
predicho = c()
for (i in (1: length(datosX))){
predicho[i] = promediosMoviles(datosX, datosY, (datosX[i]), h)
}
return(predicho)
}
yhat = (generarColumnaPrediccionesPromediosMoviles(dftest$Sunshine, dftest$RainTomorrow, 0.8))
res = sum(abs(yhat-dftest$RainTomorrow))/length(yhat)
print(paste("El porcentaje de error es ", res*100, "%") )
[1] "El porcentaje de error es 18.6868686868687 %"
LS0tCnRpdGxlOiAiVFAgMiBDbGFzaWZpY2FjaW9uIgpvdXRwdXQ6CiAgaHRtbF9ub3RlYm9vazogZGVmYXVsdAogIHBkZl9kb2N1bWVudDogZGVmYXVsdAotLS0KVFAgMiBDbGFzaWZpY2FjaW9uLgpUb21hcyBQYWxhenpvIHkgQXhlbCBGcmlkbWFuCgpDYXJnYW1vcyBsYXMgbGlicmVyaWFzCmBgYHtyfQpsaWJyYXJ5KCJnZ3Bsb3QyIikgICAgICAgICAgICAgICAgICAKbGlicmFyeSgiR0dhbGx5IikKYGBgCgpgYGB7cn0KZGYgPSByZWFkLmNzdigibGx1dmlhQXVzLmNzdiIpCmBgYApTZXRlbyBkZSBzZW1pbGxhIHJhbmRvbQpgYGB7cn0Kc2V0LnNlZWQoMzI5KQpgYGAKCkxpbXBpZXphIHkgcHJlIHByb2Nlc2FtaWVudG8gZGUgZGF0b3M6CmBgYHtyfQpkZiRSYWluVG9kYXkgPSBhcy5mYWN0b3IoZGYkUmFpblRvZGF5KSAjIHBhc28gdmFyaWFibGVzIGNhdGVnb3JpY2FzIGNvbW8gZmFjdG9yCmRmJFJhaW5Ub21vcnJvdyA9IGFzLmZhY3RvcihkZiRSYWluVG9tb3Jyb3cpCgpkZiRYIDwtIE5VTEwgIyBib3JybyBjb2x1bW5hIFggeWEgcXVlIHNvc3BlY2hhbW9zIHF1ZSBubyByZXByZXNlbnRhIG5hZGEgc2lubyBxdWUgZXMgYWxndW4gdGlwbyBkZSAiaWQiIHF1ZSBxdWVkbyBncmFiYWRvIGVuIGVsIGRhdGFmcmFtZSB5IG5vIHRpZW5lIGluZmx1ZW5jaWEgZW4gbGEgb2JzZXJ2YWNpb24uIApgYGAKCkNoZXF1ZWFtb3MgcXVlIGNhZGEgY29sdW1uYSBzZWEgZGVsIHRpcG8gY29ycmVjdG8KYGBge3J9CnN0cihkZikgCmBgYApEaXZpZG8gZGF0YXNldCBlbiB0cmFpbiB5IHRlc3QKYGBge3J9CnNhbXBsZSA8LSBzYW1wbGUoYyhUUlVFLCBGQUxTRSksIG5yb3coZGYpLCByZXBsYWNlPVRSVUUsIHByb2I9YygwLjgsMC4yKSkKZGZ0cmFpbiAgPC0gZGF0YS5mcmFtZShkZltzYW1wbGUsIF0pCmRmdGVzdCAgIDwtIGRhdGEuZnJhbWUoZGZbIXNhbXBsZSwgXSkKYGBgCgpBbmFsaXNpcyBleHBsb3JhdG9yaW8gZGUgZGF0b3MgeSB2aXN1YWxpemFjaW9uZXM6CmBgYHtyIGVjaG89VFJVRSwgZmlnLmhlaWdodD0yMCwgZmlnLndpZHRoPTIwLCBtZXNzYWdlPUZBTFNFLCB3YXJuaW5nPUZBTFNFfQpnID0gZ2dwYWlycyhkZiwgcHJvZ3Jlc3MgPSBGQUxTRSwgYmlucz0xMCkrdGhlbWVfYncoKQpnCmBgYAoKTG8gcXVlIG9ic2VydmFtb3MgZXMgcXVlIGhheSB2YXJpYWJsZXMgcXVlIGVzdGFuIG11eSBjb3JyZWxhY2lvbmFkYXMgbGluZWFsbWVudGUsIGNvbW8gbGEgZGUgdGVtcGVyYXR1cmEgbWF4aW1hIHkgdGVtcGVyYXR1cmEgYSBsYXMgM3BtLiBZIG90cmFzIHF1ZSBwYXJlY2VuIHNlciBtdXkgcG9jbyByZWxhY2lvbmFkYXMgY29tbyBsYSBodW1lZGFkIHkgbGEgcHJlc2lvbi4gCmBgYHtyfQpwcmludChwYXN0ZSgiQ29ycmVsYWNpb24gaHVtZWRhZCB5IHByZXNpb24gKGJhamEpICIgLCBjb3IoZGYkSHVtaWRpdHkzcG0sIGRmJFByZXNzdXJlM3BtKSApKQpwcmludChwYXN0ZSgiQ29ycmVsYWNpb24gbWF4VGVtcCB5IHRlbXAgYSAzcG0gKG11eSBhbHRhKSIgLCBjb3IoZGYkVGVtcDNwbSwgZGYkTWF4VGVtcCkgKSkKCmBgYAoKVGFtYmllbiBub3RhbW9zIHF1ZSBhcHJveGltYWRhbWVudGUgNC81IGRlIGxhcyBvYnNlcnZhY2lvbmVzIG5vIGxsdWV2ZSwgdGFudG8gZXNlIG1pc21vIGRpYSBjb21vIGVsIHNpZ3VpZW50ZS4KYGBge3J9Cgp0YWJsZShkZiRSYWluVG9kYXkpCnRhYmxlKGRmJFJhaW5Ub21vcnJvdykKYGBgCgpQZXJvIG5vIHNvbiBwYXJhIG5hZGEgaW5kZXBlbmRpZW50ZXMgc2kgbGx1ZXZlIGhveSB5IHNpIGxsdWV2ZSBtYcOxYW5hLiBZYSBxdWUsIHNpIHNvbG8gdHV2aWVyYSBkZSBpbmZvcm1hY2lvbiBlc3RhcyAyIGNvbHVtbmFzLCBkYWRvIHF1ZSBsbG92aW8gaG95IGxhIG51ZXZhIHByb2JhYmlsaWRhZCBkZSBxdWUgbGx1ZXZhIG1hw7FhbmEgZXMgYXByb3hpbWFkYW1lbnRlIDY1Lyg5Nys2NSkgPWFwcm94IDQwJSBtdWNobyBtYXMgcXVlIGVsIDIwJSBuYWl2ZS4gCmBgYHtyfQpkZkxsdXZpYXMgPSBkZltjKCJSYWluVG9kYXkiLCAiUmFpblRvbW9ycm93IildCnRhYmxlKGRmTGx1dmlhcykKYGBgCgpFamVyY2ljaW8gMgpgYGB7cn0KZ2dwbG90KGRmLCBhZXMoeD1TdW5zaGluZSwgeT1IdW1pZGl0eTNwbSwgY29sb3I9UmFpblRvbW9ycm93KSkgKwogIGdlb21fcG9pbnQoKSAKYGBgCk9ic2VydmFtb3MgY29zYXMgbXV5IHJlbGV2YW50ZXMsIGxvcyBkaWFzIHF1ZSB2YSBhIGxsb3ZlciBtYcOxYW5hLCB0aWVuZW4gbWF5b3IgaHVtZWRhZCB5IG1lbm9zIHNvbCwgeSBsb3MgZGlhcyBxdWUgbm8gbGxvdmVyYSBtYcOxYW5hIHRpZW5lbiB0b2RvcyBtdWNobyBtYXMgc29sIHkgbGEgaHVtZWRhZCBlbiBwcm9tZWRpbyBlcyBtYXMgYmFqYS4gQSBzdSB2ZXogaGF5IHZhcmlvcyBkaWFzLCBlbiBzdSBtYXlvcmlhIGRpYXMgcXVlIGxsb3ZlcmEgbWHDsWFuYSwgY3V5byBuaXZlbCBkZSBzb2wgZXMgMCwgbG8gY3VhbCBnZW5lcmEgZXNhIGNvbHVtbmEgZW4gZWwgbGFkbyBpenF1aWVyZG8uIApUYW1iaWVuIHZlbW9zIHF1ZSBzaSBiaWVuIGhheSBtdWNob3MgZGlhcyBxdWUgdGllbmVuIG11Y2hvIHNvbCB5IHBvY2EgaHVtZWRhZCAobG9zIGRpYXMgcXVlIG5vIGxsb3ZlcmEgbWHDsWFuYSksIG5vIHZlbW9zIGNhc2kgbmluZ3VuYSBvYnNlcnZhY2lvbiBjb24gcG9jbyBzb2wgeSBwb2NhIGh1bWVkYWQsIGxvIGN1YWwgbm9zIHBvZHJpYSBoYWJsYXIgZGUgY2llcnRhIHJlbGFjaW9uIGh1bWVkYWQgLSBzb2wuIApOb3NvdHJvcyBwZW5zYW1vcyBxdWUgbGEgaHVtZWRhZCB0aWVuZSBtYXlvciBjYXBhY2lkYWQgcHJlZGljdGl2YSAoc2kgc2UgdG9tYXJhIHVuYSBzb2xhIHkgbm8gZW4gY29uanVudG8gYW1iYXMpLCB5YSBxdWUgbGEgdmFyaWFibGUgc3Vuc2hpbmUgZXN0YSBtdWNobyBtYXMgZGlzcGVyc2EgcGFyYSBsb3MgZGlhcyBxdWUgTm8gbGxvdmlvIGFsIGRpYSBzaWd1aWVudGUsIGxvIGFuYWxpemFyZW1vcyBlbiBkb3MgZ3JhZmljb3MgZGUgZGVuc2lkYWQuCgpgYGB7cn0KcHM8LWdncGxvdChkZiwgYWVzKHg9U3Vuc2hpbmUsIGZpbGw9UmFpblRvbW9ycm93KSkgKwogIGdlb21fZGVuc2l0eShhbHBoYT0wLjQpICsgbGFicyh4PSAiTml2ZWwgZGUgcmFkaWFjaW9uIHNvbGFyIChTdW5zaGluZSkiLAogICAgICAgc3VidGl0bGU9IkdyYWZpY28gZGVuc2lkYWQgcmFkaWFjaW9uIHNvbGFyIikgKyBnZW9tX3ZsaW5lKHhpbnRlcmNlcHQ9Ny44LCBzaXplPTAuNSwgY29sb3I9InJlZCIpCnBzCmBgYApgYGB7cn0KcGg8LWdncGxvdChkZiwgYWVzKHg9SHVtaWRpdHkzcG0sIGZpbGw9UmFpblRvbW9ycm93KSkgKwogIGdlb21fZGVuc2l0eShhbHBoYT0wLjQpICsgbGFicyh4PSAiTml2ZWwgZGUgaHVtZWRhZCAoSHVtaWRpdHkzcG0pIiwKICAgICAgIHN1YnRpdGxlPSJHcmFmaWNvIGRlbnNpZGFkIGh1bWVkYWQiKSArIGdlb21fdmxpbmUoeGludGVyY2VwdD02Miwgc2l6ZT0wLjUsIGNvbG9yPSJyZWQiKQoKcGgKYGBgCkRlc3B1ZXMgZGUgdmVyIGVzdG9zIDIgZ3JhZmljb3Mgbm90YW1vcyBxdWUgbm8gZXMgZmFjaWwgZGFyIHVuIHB1bnRvIGRlIGNvcnRlIHBhcmEgZGlmZXJlbmNpYXIgbGFzIDIgY2xhc2VzIHNvbGFtZW50ZSB0b21hbmRvIHVuYSB2YXJpYWJsZS4gWWEgcXVlIHNpIHRvbWFzZW1vcyBhcHJveGltYWRhbWVudGUgNy44IGRlIHB1bnRvIGRlIGNvcnRlIHBhcmEgbGEgcmFkaWFjaW9uIHNvbGFyIG8gNjIgcGFyYSBsYSBodW1lZGFkIGNvbW8gcHVudG8gZGUgY29ydGUsIGRlIHRvZGFzIGZvcm1hcyB0ZW5kcmlhcyBiYXN0YW50ZSBlcnJvciB5YSBxdWUgbGFzIGNsYXNlcyBzZSBzb2xhcGFuIG11Y2hvIG1pcmFuZG9sYXMgdW5pZGltZW5zaW9uYWxtZW50ZS4gVG9tYW1vcyBlc3RvcyB2YWxvcmVzIGRlIHJlZmVyZW5jaWEgY29tbyBwYXJhIGRlY2lyIHF1ZSBuaW5ndW4gY29ydGUgZXMgYnVlbm8sIGVzdG9zIG5pIHNpcXVpZXJhIHRpZW5lbiBlbiBjdWVudGEgbGEgZGlmZXJlbmNpYSBkZSBwcm9wb3JjaW9uIGRlIGNsYXNlcy4gRW4gZGVmaW5pdGl2YSBubyBub3MgY2FzYW1vcyBjb24gbmluZ3VuYSB2YXJpYWJsZS4gCgpFamVyY2ljaW8gMy4KTG9zIGJveHBsb3Qgbm8gc29uIGJ1ZW5vcyBncmFmaWNvcy4gUHVlZGVuIG9jdWx0YXIgZGVtYXNpYWRhIGluZm9ybWFjaW9uIGN1YW5kbyBob3kgZW4gZGlhIHRlbmVtb3MgbGEgY2FwYWNpZGFkIGRlIHByb2Nlc2FybGEuIAoKYGBge3J9CnAyPC1nZ3Bsb3QoZGYsIGFlcyh5PUh1bWlkaXR5M3BtLCB4PVJhaW5Ub21vcnJvdywgZmlsbD1SYWluVG9tb3Jyb3cpKSArCiAgZ2VvbV9ib3hwbG90KCkgKyBsYWJzKHg9ICJOaXZlbCBkZSBodW1lZGFkIChIdW1pZGl0eTNwbSkiLAogICAgICAgc3VidGl0bGU9IkJveHBsb3RzIGh1bWVkYWQgc2VndW4gbGx1dmlhIG1hw7FhbmEiKQpwMgpgYGAKYGBge3J9CnAzPC1nZ3Bsb3QoZGYsIGFlcyh5PVN1bnNoaW5lLCB4PVJhaW5Ub21vcnJvdywgZmlsbD1SYWluVG9tb3Jyb3cpKSArCiAgZ2VvbV9ib3hwbG90KCkgKyBsYWJzKHg9ICJOaXZlbCBkZSByYWRpYWNpb24gc29sYXIgKFN1bnNoaW5lKSIsCiAgICAgICBzdWJ0aXRsZT0iQm94cGxvdHMgcmFkaWFjaW9uIHNvbGFyIHNlZ3VuIGxsdXZpYSBtYcOxYW5hIikKcDMKYGBgCkNvbW8gdmVuIGxhcyBtZWRpYXMgZGlmaWVyZW4gcGFyYSBhbWJvcyBjYXNvcywgcGVybyBlc2EgaW5mb3JtYWNpb24geWEgbGEgaGFiaWFtb3MgdmlzdG8gKGFkZW1hcyBkZSBtdWNoYXMgb3RyYXMgY29zYXMgcXVlIGFjYSBubykgZW4gbG9zIGRlbnNpdHkgcGxvdHMuIE5vIGhheSBvdXRsaWVycyAvIHZhbG9yZXMgYXRpcGljb3MuIAoKRWplcmNpY2lvIDQKUGFyYSBoYWNlciBsYXMgdmVudGFuYXMgbW92aWxlcyB2b3kgYSBwcmltZXJvIHRyYW5zZm9ybWFyIGVsIGRhdGFzZXQgZW4gMSB5IDAgYSBsYXMgY2F0ZWdvcmljYXMsIHBhcmEgcG9kZXIgbHVlZ28gdG9tYXIgcHJvbWVkaW9zLgpgYGB7cn0KZGZ0cmFpbiRSYWluVG9kYXkgPSBpZmVsc2UoZGZ0cmFpbiRSYWluVG9kYXk9PSJZZXMiLDEsMCkKZGZ0cmFpbiRSYWluVG9tb3Jyb3cgPSBpZmVsc2UoZGZ0cmFpbiRSYWluVG9tb3Jyb3c9PSJZZXMiLDEsMCkKZGZ0ZXN0JFJhaW5Ub2RheSA9IGlmZWxzZShkZnRlc3QkUmFpblRvZGF5PT0iWWVzIiwxLDApCmRmdGVzdCRSYWluVG9tb3Jyb3cgPSBpZmVsc2UoZGZ0ZXN0JFJhaW5Ub21vcnJvdz09IlllcyIsMSwwKQpgYGAKCmBgYHtyfQpwcm9tZWRpb3NNb3ZpbGVzID0gZnVuY3Rpb24oZGF0b3NYLCBkYXRvc1ksIHZhbG9yLCBoKXsKICBkZmUgPSBkYXRhLmZyYW1lKGRhdG9zWCxkYXRvc1kpCiAgZGYyID0gZGZlW2RmZSRkYXRvc1ggPj0gdmFsb3ItaCAmIGRmZSRkYXRvc1ggPD0gdmFsb3IraCAsXQogIHdoaWxlKG5yb3coZGYyKT09MCl7CiAgICBoID0gaCoyICMgU2kgbm8gbWUgYWdhcnJhIGEgbmFkaWUgcGFyYSBwcm9tZWRpYXIsIGF1bWVudGFtZSBsYSB2ZW50YW5pdGEuIE5vIHN1ZWxlIHBhc2FyLCBzb2xvIG91dGxpZXJzIG8gZXh0cmVtb3MKICAgIGRmMiA9IGRmZVtkZmUkZGF0b3NYID49IHZhbG9yLWggJiBkZmUkZGF0b3NYIDw9IHZhbG9yK2ggLF0KICB9CiAgaWYobWVhbihkZjIkZGF0b3NZKT4xLzIgKXsKICAgIHJldHVybigxKQogIH0KICByZXR1cm4oMCkKfQpgYGAKYGBge3J9CnByb21lZGlvc01vdmlsZXMoZGZ0cmFpbiRTdW5zaGluZSwgZGZ0cmFpbiRSYWluVG9tb3Jyb3csIDgsIDAuMSkKYGBgCmBgYHtyfQpwcm9tZWRpb3NNb3ZpbGVzKGRmdHJhaW4kSHVtaWRpdHkzcG0sIGRmdHJhaW4kUmFpblRvbW9ycm93LCA4NSwgMikKYGBgCgpFamVyY2ljaW8gNSAKTm9zIGNyZWFtb3MgdW5hIGZ1bmNpb24gcXVlIG5vcyBnZW5lcmUgdG9kbyBlbCB2ZWN0b3IgZGUgcHJlZGljY2lvbmVzIFloYXQuCgpWYW1vcyBhIGV2YWx1YXJsbyBjb24gZWwgbWV0b2RvIGRlIHZhbGlkYWNpb24gY3J1emFkYSBkZSBMT08gKGRlamFyIHVubyBhZnVlcmEgcGFyYSBlbnRyZW5hciB5IGV2YWx1YXJsbyBjb24gZXNlKS4KCmBgYHtyfQpsZWF2ZU9uZU91dCA9IGZ1bmN0aW9uKGRhdG9zWCwgZGF0b3NZLCBoKXsKICBlcnJvciA9IDAKICBmb3IgKGkgaW4gKDE6IGxlbmd0aChkYXRvc1gpKSl7CiAgICBwcmVkaWNob0kgPSBwcm9tZWRpb3NNb3ZpbGVzKGRhdG9zWFstaV0sIGRhdG9zWVstaV0sIGRhdG9zWFtpXSwgaCkKICAgIGVycm9yID0gZXJyb3IgKyBhYnMocHJlZGljaG9JIC0gZGF0b3NZW2ldKQogIH0KICByZXR1cm4oZXJyb3IpCn0KYGBgCgpgYGB7cn0KaFBvc2libGVIdW1lZGFkID0gc2VxKDEsIDMwLCAwLjUgKQpoUG9zaWJsZVN1bnNoaW5lID0gc2VxKDAuNSwgMTAsIDAuMSApCmBgYAoKCmBgYHtyfQplcnJvcmVzaEh1bWVkYWQgPSBjKCkKaHVtID0gZGZ0cmFpbiRIdW1pZGl0eTNwbQpsbHV2ID0gZGZ0cmFpbiRSYWluVG9tb3Jyb3cKZm9yIChpIGluICgxOiBsZW5ndGgoaFBvc2libGVIdW1lZGFkKSkpewogICAgZXJyb3Jlc2hIdW1lZGFkW2ldID0gbGVhdmVPbmVPdXQoaHVtLCBsbHV2LCBoUG9zaWJsZUh1bWVkYWRbaV0pCn0KYGBgCgpgYGB7cn0KZXJyb3Jlc2hTdW5zaGluZSA9IGMoKQpodW0gPSBkZnRyYWluJEh1bWlkaXR5M3BtCmxsdXYgPSBkZnRyYWluJFJhaW5Ub21vcnJvdwpzdW4gPSBkZnRyYWluJFN1bnNoaW5lCmZvciAoaSBpbiAoMTogbGVuZ3RoKGhQb3NpYmxlU3Vuc2hpbmUpKSl7CiAgICBlcnJvcmVzaFN1bnNoaW5lW2ldID0gbGVhdmVPbmVPdXQoc3VuLCBsbHV2LCBoUG9zaWJsZVN1bnNoaW5lW2ldKQp9CmBgYApgYGB7cn0KcGxvdChoUG9zaWJsZUh1bWVkYWQgLCBlcnJvcmVzaEh1bWVkYWQsIHR5cGUgPSAibCIpCmBgYAoKRGUgYWNhIHZlbW9zIHF1ZSBsYSB2ZW50YW5hIG9wdGltYSBwYXJhIGh1bWVkYWQgZXMgMS4KCmBgYHtyfQpwbG90KGhQb3NpYmxlU3Vuc2hpbmUgLCBlcnJvcmVzaFN1bnNoaW5lLCB0eXBlID0gImwiKQpgYGAKTWllbnRyYXMgcXVlIGxhIHZlbnRhbmEgb3B0aW1hIHBhcmEgbGEgcmFkaWFjaW9uIHNvbGFyIGVzIGRlIDAuOAoKRWplcmNpY2lvIDYKYWhvcmEgbm9zIGZpamFtb3MgZWwgZXJyb3IgY29uIGVsIGRhdGFzZXQgZGUgdGVzdCBxdWUgc2VwYXJhbW9zIGFsIHByaW5jaXBpbyBwYXJhIGVzdGEgdmVudGFuYQpgYGB7cn0KZ2VuZXJhckNvbHVtbmFQcmVkaWNjaW9uZXNQcm9tZWRpb3NNb3ZpbGVzID0gZnVuY3Rpb24oZGF0b3NYLCBkYXRvc1kgLCBoKXsKICBwcmVkaWNobyA9IGMoKQogIGZvciAoaSBpbiAoMTogbGVuZ3RoKGRhdG9zWCkpKXsKICAgIHByZWRpY2hvW2ldID0gcHJvbWVkaW9zTW92aWxlcyhkYXRvc1gsIGRhdG9zWSwgKGRhdG9zWFtpXSksIGgpCiAgfQogIHJldHVybihwcmVkaWNobykKfQpgYGAKYGBge3J9CnloYXQgPSAoZ2VuZXJhckNvbHVtbmFQcmVkaWNjaW9uZXNQcm9tZWRpb3NNb3ZpbGVzKGRmdGVzdCRTdW5zaGluZSwgZGZ0ZXN0JFJhaW5Ub21vcnJvdywgMC44KSkKYGBgCmBgYHtyfQpyZXMgPSBzdW0oYWJzKHloYXQtZGZ0ZXN0JFJhaW5Ub21vcnJvdykpL2xlbmd0aCh5aGF0KQpwcmludChwYXN0ZSgiRWwgcG9yY2VudGFqZSBkZSBlcnJvciBlcyAiLCByZXMqMTAwLCAiJSIpICkKYGBgCgo=